WindFluxInit Subroutine

SUBROUTINE WindFluxInit &
!
( ini, mask, dtMeteo, tstart, dem, dem_loaded)

IMPLICIT NONE

TYPE (IniList), INTENT(IN) :: ini
TYPE (grid_integer), INTENT(IN) :: mask  !!defines interpolation extent
INTEGER (KIND = short), INTENT(IN) :: dtMeteo !! deltat of meteo data reading
TYPE (DateTime),     INTENT(IN) :: tstart !!initial time
TYPE(grid_real),    INTENT(in)  :: dem !!digital elevation model to be used to modify interpolated data
LOGICAL , INTENT (IN) :: dem_loaded !! true if dem has been loaded


!local declarations
CHARACTER (LEN = 1000) :: filename, filenameWD
TYPE (grid_real) :: meteoTemp


!-------------------------end of declarations----------------------------------
!set initial time
  timeNew = tstart  

  !initialize grid
  CALL NewGrid (windSpeed, mask, 0.)
  CALL NewGrid (grid_devst, mask, 0.)
  
  !set deltat
  dtWindSpeed = IniReadInt ('dt', ini, section = 'wind-speed')

  !check dt is multiple of dtMeteo
  IF (.NOT.(MOD(dtWindSpeed,dtMeteo) == 0)) THEN
           CALL Catch ('error', 'WindFlux',   &
		                 'dt is not multiple of dtMeteo')
  END IF
  
  !set valid threshold
  IF (KeyIsPresent ('valid-threshold', ini, &
                         section = 'wind-speed') ) THEN
      valid_prcn = IniReadReal ('valid-threshold', ini, &
                                        section = 'wind-speed')
  ELSE ! set to default = 1.0 
      CALL Catch ('info', 'WindFlux',   &
		                 'valid-threshold not defined, set to 1')
      valid_prcn = 1.
  END IF

  !set cell-size
  cellsizeInterpolation = mask % cellsize
  
  !interpolation-assignment method
  IF (KeyIsPresent('interpolation-assignment', ini, &
                     section = 'wind-speed') ) THEN
      interpolationMethod_assignment = IniReadInt ('interpolation-assignment', ini, &
                                    section = 'wind-speed')
  ELSE
    CALL Catch ('error', 'WindFlux',   &
                'interpolation-assignment missing in meteo configuration file')
  END IF   

  !set interpolation method
  IF (interpolationMethod_assignment == 1) THEN

      interpolationMethod = IniReadInt ('interpolation', ini, &
                                    section = 'wind-speed')
      
	    CALL SetSpecificProperties (interpolationMethod, ini, mask)
  
  ELSE !read map
      interpolationMethod = -1
      CALL GridByIni (ini, interpolationMethod_map, 'wind-speed', &
                     'interpolation')
      !scan for interpolation methods included
      interpolationMethod_vector = 0
      DO i = 1, interpolationMethod_map % idim
        DO j = 1, interpolationMethod_map % jdim
            IF (interpolationMethod_map % mat(i,j) /= &
                interpolationMethod_map % nodata) THEN
                interpolationMethod_vector &
                   (interpolationMethod_map % mat (i,j)) = &
                    interpolationMethod_map % mat (i,j)
            END IF
        END DO
      END DO

      DO i = 1, 5
         CALL SetSpecificProperties (interpolationMethod_vector (i), ini, mask)
      END DO

  END IF

  !scale factor and offset
	IF (KeyIsPresent('offset', ini, section = 'wind-speed')) THEN	
	   offset_value = IniReadReal ('offset', ini, section = 'wind-speed')
	ELSE
	    offset_value = MISSING_DEF_REAL
	END IF
	
	IF (KeyIsPresent('scale-factor', ini, section = 'wind-speed')) THEN	
	   scale_factor = IniReadReal ('scale-factor', ini, section = 'wind-speed')
	ELSE
	    scale_factor = MISSING_DEF_REAL
	END IF			

   !set power_idw
   IF (KeyIsPresent('idw-power', ini, section = 'wind-speed')) THEN	
	   idw_power = IniReadReal ('idw-power', ini, section = 'wind-speed')
   ELSE !set default value
	    idw_power = 2.
   END IF	
 
	!file
    filename = IniReadString ('file', ini, section = 'wind-speed')
    IF (interpolationMethod_assignment == 1 .AND. &
        interpolationMethod == 0) THEN !data are stored in net-cdf file
       !store net-cdf filename
       fileGrid = filename
       IF ( KeyIsPresent ('variable', ini, section = 'wind-speed') ) THEN
          windSpeed % var_name = IniReadString ('variable', &
                                             ini, section = 'wind-speed')
            !read grid and store as temporary file
             CALL NewGrid (layer = meteoTemp, fileName = TRIM(fileGrid), &
                         fileFormat = NET_CDF, &
                         variable = TRIM(windSpeed % var_name) )
       ELSE IF  (KeyIsPresent ('standard_name', ini, &
                              section = 'wind-speed') ) THEN
          windSpeed % standard_name = IniReadString ('standard_name', &
                                              ini, section = 'wind-speed')
       ELSE
          CALL Catch ('error', 'WindFlux',   &
		       'standard_name or variable missing in section wind-speed' )
       END IF      

       !set cellsize to zero. Optimal cellsize is automatically computed
       cellsizeInterpolation = 0.
            
    ELSE !open file containing local measurements
       fileunit = GetUnit ()
	     OPEN(fileunit,file = filename(1:LEN_TRIM(filename)),status='old')
    END IF
    
    elevationDrift = 0
    micrometInUse = .FALSE.
    !set elevationfridt option
    IF ( interpolationMethod_assignment == 1 ) THEN
        
        IF ( interpolationMethod == 4 ) THEN  !micromet
           elevationDrift = 1
           micrometInUse = .TRUE.
        END IF
        
        IF ( interpolationMethod == 5 ) THEN  !gonzalez
            elevationDrift = 1
        END IF
            
    ELSE !interpolationMethod_assignment = 2
        !parse interpolationMethod_map
        DO i = 1, interpolationMethod_map % idim
            DO j = 1, interpolationMethod_map % jdim
                IF ( interpolationMethod_map % mat (i,j) == 4 ) THEN ! micromet
                     elevationDrift = 1
                     micrometInUse = .TRUE.
                END IF
                IF ( interpolationMethod_map % mat (i,j) == 5 ) THEN !gonzalez
                    elevationDrift = 1
                END IF
            END DO
        END DO 
    END IF
      
    IF (elevationDrift == 1 ) THEN
      
          !check if dem have been loaded by domain properties
          IF ( .NOT. dem_loaded) THEN
               CALL Catch ('error', 'WindFlux',   &
                            'dem for elevation drift was not loaded ')
          END IF
          
          
          !load wind direction metadata
         IF (KeyIsPresent ('wind-direction-file', ini, &
                         section = 'wind-speed') ) THEN
               filenameWD = IniReadString ('wind-direction-file', &
                                       ini, section = 'wind-speed')
        ELSE
             CALL Catch ('error', 'WindFlux',   &
		                 'missing wind speed direction file &
                          in section wind-speed')                
        
        END IF
          
        fileunitWD = GetUnit ()
        
        OPEN (fileunitWD, file = filenameWD(1:LEN_TRIM(filenameWD)), status='old')

        CALL ReadMetadata (fileunitWD, anemometersWD)
          
    END IF
    
    IF ( micrometInUse ) THEN
       
       !load micromet parameters values
       IF (KeyIsPresent ('micromet-length-scale', ini, &
                         section = 'wind-speed') ) THEN
           micrometLengthScale = IniReadReal ('micromet-length-scale', ini, &
                                        section = 'wind-speed')
       ELSE ! set to default = 5000
            CALL Catch ('info', 'WindFlux',   &
		                 'micrometLengthScale set to 5000 m')
            micrometLengthScale = 5000.
       END IF
       
       IF (KeyIsPresent ('micromet-slopewt', ini, &
                         section = 'wind-speed') ) THEN
           micrometSlopeWT = IniReadReal ('micromet-slopewt', ini, &
                                        section = 'wind-speed')
       ELSE ! set to default = 0.5
            CALL Catch ('info', 'WindFlux',   &
		                 'micrometSlopeWT set to 0.5')
            micrometSlopeWT = 0.5
       END IF
       
       IF (KeyIsPresent ('micromet-curvewt', ini, &
                         section = 'wind-speed') ) THEN
           micrometCurvatureWT = IniReadReal ('micromet-curvewt', ini, &
                                        section = 'wind-speed')
       ELSE ! set to default = 0.5
            CALL Catch ('info', 'WindFlux',   &
		                 'micrometCurvatureWT set to 0.5')
            micrometCurvatureWT = 0.5
       END IF
       
       !Compute the curvature
       CALL CurvatureMicroMet (dem, micrometLengthScale, curvature)

       !Compute the terrain slope [rad]
       CALL DeriveSlope (dem, slope)
  
       ! Compute the terrain slope azimuth [rad]
       CALL DeriveAspect (dem, slope_az)
       
        
    END IF  

  !grid exporting settings
  IF (KeyIsPresent ('export', ini, section = 'wind-speed')  )  THEN
     export = IniReadInt ('export', ini, section = 'wind-speed')
  ELSE
     export = 0
  END IF 

  IF (export == 1) THEN
     
     
     !set export path name
     IF (KeyIsPresent ('export-path', ini, section = 'wind-speed')  )  THEN
         export_path = IniReadString ('export-path', ini, section = 'wind-speed')
         !check if path ends with / or \
         IF (  export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '\' .AND. &
               export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '/' ) THEN
             export_path (LEN_TRIM (export_path)+1:LEN_TRIM (export_path)+1) = '\'
         END IF
       
         IF (INDEX (export_path,'.') == 1) THEN !detected relative path, remove '.'
            export_path = export_path (2:LEN_TRIM(export_path))
            !build absolute path
            export_path = TRIM(CurrentDir() ) // TRIM(export_path)
         END IF

         !check OS convention
        IF (GetOS () == WIN32) THEN
          export_path = ReplaceChar (export_path,'/','\')
        ELSE
          export_path = ReplaceChar (export_path,'\','/')
        END IF
  
         
         !check folder exists
         IF ( .NOT. DirExists (TRIM (export_path) ) ) THEN
             CALL Catch ('info', 'WindFlux',   &
                  'creating directory:  ', argument = TRIM(export_path))
             CALL DirNew (export_path)
         END IF
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-path')
     END IF 
     
     !starting time
     IF (KeyIsPresent ('export-start', ini, section = 'wind-speed')  )  THEN
         timeString = IniReadString ('export-start', ini, 'wind-speed')
         export_start = timeString
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-start')
     END IF 
     
     !initialize timeNewExport
     timeNewExport = export_start
     
     !ending time
     IF (KeyIsPresent ('export-stop', ini, section = 'wind-speed')  )  THEN
         timeString = IniReadString ('export-stop', ini, 'wind-speed')
         export_stop = timeString
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-start')
     END IF 
     
     !export dt
     IF (KeyIsPresent ('export-dt', ini, section = 'wind-speed')  )  THEN
         export_dt = IniReadInt ('export-dt', ini, section = 'wind-speed')
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-dt')
     END IF 

     !coordinate reference system
     IF (KeyIsPresent ('export-epsg', ini, section = 'wind-speed')  )  THEN
         export_epsg = IniReadInt ('export-epsg', ini, section = 'wind-speed')
         
         exportGridMapping = DecodeEPSG (export_epsg)
         IF (exportGridMapping == windSpeed % grid_mapping) THEN
            needConversion = .FALSE.
            !initialize grid for exporting with CRS as meteo variable
            CALL NewGrid (exportedGrid, windSpeed)
            CALL NewGrid (exportedGridVar, windSpeed)
         ELSE
            needConversion = .TRUE.
            !initialize grid for converting with required CRS
            exportedGrid % grid_mapping = DecodeEPSG (export_epsg)
            exportedGridVar % grid_mapping = DecodeEPSG (export_epsg)
         END IF
     ELSE
         CALL Catch ('info', 'WindFlux',   &
                  'export-epsg not defined, use default')
         needConversion = .FALSE.
         !initialize grid for exporting with CRS 
         CALL NewGrid (exportedGrid, windSpeed)
         CALL NewGrid (exportedGridVar, windSpeed)
     END IF 

     !export file format 
     IF (KeyIsPresent ('export-format', ini, section = 'wind-speed')  )  THEN
         export_format = IniReadInt ('export-format', ini, section = 'wind-speed')
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-format')
     END IF 

     IF (export_format == 3) THEN
       !grid map  
       CALL SetLongName ( 'wind_speed', exportedGrid)
       CALL SetStandardName ( 'wind_speed', exportedGrid)
       CALL SetUnits ('m/s', exportedGrid) 
       !if file exists, remove it
       export_file = TRIM(export_path) //  'wind_speed.nc'
       IF ( FileExists (export_file) ) THEN
          CALL FileDelete (export_file)
       END IF
       
       !variance of kriging 
       CALL SetLongName ( 'wind_speed_variance', exportedGridVar)
       CALL SetStandardName ( 'wind_speed_variance', exportedGridVar)
       CALL SetUnits ('m/s', exportedGridVar) !this unit is for exporting grid
       !if file exists, remove it
       export_file_var = TRIM(export_path) //  'wind_speed_variance.nc'
       IF ( FileExists (export_file_var) ) THEN
          CALL FileDelete (export_file_var)
       END IF
       
     END IF

  END IF
	
	!complete initialization

  IF (interpolationMethod == 0) THEN
		!Get the dt of imported  field. Assume dt is regular	
		dtGrid = GetDtGrid (filename = fileGrid, checkRegular = .TRUE.)
		!check dt is multiple of dtGrid
        
        IF (.NOT.(MOD(dtWindSpeed,dtGrid) == 0)) THEN
            CALL Catch ('error', 'WindFlux',   &
                'dt wind speed is not multiple of dt of input grid')
        END IF
   ELSE
        !populate  metadata
        CALL ReadMetadata (fileunit, anemometersWS)

        !check spatial reference system
        IF ( .NOT. anemometersWS % mapping == mask % grid_mapping)  THEN
            CALL Catch ('info', 'WindFlux',   &
		              'converting coordinate of stations')
            !convert stations' coordinate
            point1 % system = anemometersWS % mapping
            point2 % system = mask % grid_mapping
            anemometersWS % mapping = mask % grid_mapping
            DO i = 1, anemometersWS % countObs
              point1 % easting = anemometersWS % obs (i) % xyz % easting
              point1 % northing = anemometersWS % obs (i) % xyz % northing
              point1 % elevation = anemometersWS % obs (i) % z
              CALL Convert (point1, point2)
              anemometersWS % obs (i) % xyz % easting = point2 % easting
              anemometersWS % obs (i) % xyz % northing = point2 % northing 
            END DO
        END IF
        
       
  END IF
   
RETURN
END SUBROUTINE WindFluxInit